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A Kalman recursive algorithm for estimating time from an ensemble of atomic clocks has been developed. 
The algorithm allows for the addition or deletion of clocks at any time, and provides automatic error detec- 
tion and correction. The observations consist of time differences between clocks and may be taken at un- 
equally spaced time points. Maximum likelihood estimates of the unknown parameters are obtained with 
confidence intervals, as well as hypothesig tests to determine whether the estimated parameters are significantly 
different from zero. The program is operational on the National Bureau of Standards' Time and Frequency 
Division's PDP 11/70. 
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1. Introduction 

Cesium beam atomic clocks have an accuracy of a few 
parts in 10 14 over a period of a day; however, they are 
not deterministic and undergo stochastic variations in both 
time and frequency. It was shown by Tryon and Jones 
[1] that the actual frequency of a clock behaves as a 
random walk which, over a time period of a day, has a 
standard deviation of less than four nanoseconds per day. 
The effect of this random walk in frequency on time 
measurements is that the individual frequency steps are 
summed, producing an integrated random walk. 

In addition, frequency has a white noise component 
which integrates into an independent random walk in time 
with a standard deviation of up to 15 nanoseconds over 
a time period of a day. Therefore, the deviation of clock 
time from true time behaves as a random walk plus the 
sum of a random walk. The possibility of a frequency drift 
also exists. 

When estimating time from an ensemble of clocks, the 
only observations possible are measurements of time dif- 
ferences between clocks. In addition, observations may 
occur at unequally spaced time points, and various types 
of errors are possible. The most common errors are read 
errors, where a single measurement is incorrect but subse- 
quent measurements are the same as if the read error did 
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not occur; time steps, where a jump in time occurs and 
the clock remains at the new position; and frequency steps, 
where the frequency takes a step and remains at the new 
position. This last error appears as a change in the rate 
of gain of the clock. 

This paper describes a state space algorithm for 
estimating time from an ensemble of cesium beam atomic 
clocks with unequally spaced observations subject to 
various errors. Clocks can be added to or deleted from 
the ensemble. Using iterative calculations over a recent 
history of the ensemble, maximum likelihood estimates 
of the unknown variances can be obtained by nonlinear 
optimization. Confidence intervals on the estimated 
parameters and test of hypotheses, such as whether 
parameters are significantly different from zero, can also 
be obtained from the likelihood function. 

2. Mathematical Model 

The model for a single clock in state space form is 



(2.1) 



where t = 1, 2, 3, .■ . . is an index of the observation 
points, x{t ) is the difference between the clock's time and 
"true" time, and yit) is the difference between the clock's 
frequency and the fundamental resonance of the cesium 
atom which defines the second, 9,192,631,770 Hz, and 
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is expressed in units of nanoseconds per day. wit) 
represents a possible frequency drift in units of 
nanoseconds/day 2 , and 6it) is the time interval between 
(-1 and t in days. eU), r](t) and a{t) are random variables 
with zero mean, uncorrelated with each other, and un- 
correlated in time. These are the input to the random 
walks, and for small d(t), their variances are proportional 
to the length of the time intervals between observations, 
d(r), 



Var {£<()} = d(() o £ * 
Var fa(t)} = 6(t) oj 
Var {ait)} = 6{t) o ff 2 



(2.2) 



If a a a = 0, the drift, w{t), does not change with time. 
The state equation for an ensemble of m clocks is ob- 
tained by concatenating the three state elements of each 
clock into a column vector of length 3m. The state tran- 
sition matrix, <b(t), is a 3m by 3m matrix consisting of 
3 by 3 blocks on the main diagonal corresponding to the 
state transition matrix in eq (1), and zero blocks off the 
main diagonal. This equation of state for m clocks will 
be written. 



Xit) =Mt) Xit-l) + Uit). 



(2.3) 



The random input vector, U(t), is of length 3m and con- 
sists of the e(t), r\{t), and ait) for each clock. The 
covariance matrix of Uit), 6it) Q, is a 3m by 3m diagonal 
matrix with diagonal elements consisting of the o £ 2 , o^ 2 
and o a 2 for each clock. While reasonable guesses are 
available for the o £ 2 and o 2 , it is necessary to obtain bet- 
ter estimates of these variances from data. The values are 
characteristics of each clock and will be different for each 
clock. It is also necessary to determine whether the o a 2 
are significantly different from zero. If the a a 2 are not 
significantly different from zero it would indicate that any 
drifts that exist are deterministic rather than random 
walks. In this case the wit) will be constant with respect 
to time and these can be tested to determine if the drifts 
are significantly different from zero. 

When observations are taken, one of the clocks is used 
as a reference clock, and the time differences between this 
clock and each of the other clocks are recorded. This gives 
an observation vector Z(t) , of length m-1, and an obser- 
vation equation of the form 



zft) = Humt) + vu), 



(2.4) 



where Hit) is a matrix indicating which clock differences 
are observed. If the first clock is the reference clock, Hit) 



is of the form 

100-100 000 
100 000-100 
100 000 000-1 



Hit) = 



(2.5) 



and has dimension (m-1 ) by 3m. Hit) can be time depen- 
dent since if observations are missing or deleted because 
of error detection, rows of Hit) are eliminated. V(t) is a 
vector of observational errors. These errors are very small 
when reading clock differences. If the data are truncated 
to the nearest nanosecond, the variances of the elements 
Vit) would be 



r = 1/12 (nanoseconds) 2 . 



(2.6) 



When observations are taken with higher precision than 
the nearest nanosecond, r can be reduced and its value 
determined from the instrumentation. The observational 
error variance matrix, R, will be a diagonal matrix and 
diagonal elements r, assumed to be equal and known. 
While these numbers are very small, the inclusion of R 
in the recursion does add numerical stability. 



3. The Data 

In a companion paper [1] Tryon and Jones used 333 
daily observations from seven clocks starting February 16, 
1979 to estimate clock parameters. These data were 
chosen since they consisted of a fairly long record with 
few errors. Some preprocessing was necessary to remove 
several outliers giving a set of data that could be con- 
sidered error-free. These calculations were carried out on 
a GDC Cyber 170/750 computer. 

This paper reports the new statistical procedures and 
techniques that were incorporated when the algorithm was 
rewritten to run on the National Bureau of Standards 
Time and Frequency Division's PDP-11/70. The 
algorithm was initialized on March 31, 1981 with the first 
set of data collected on April 1, 1981. At the beginning 
of each month the daily observations from the previous 
month are passed to the algorithm without preprocessing. 
While the data are nominally collected at the same time 
each day, for various reasons, the time of the observa- 
tions are sometimes off by several hours giving unequally 
spaced data. Individual clocks may have read errors, time 
steps, and frequency steps. In addition, clocks may be 
added or deleted from the ensemble. 
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4. The Kalman Recursion With 
Error Detection 

Let X(r|s), * < t be the best estimate of the state vector 
at time t given obervations up to time s, and let Pit \s) be 
its eovariance matrix. The recursion begins by specifying 
X(0|0) and P(0|0). First a one step prediction of the state 
vector is calculated [3] 



/<(+!) = A { p + c, 



where 



Ai'- [0 



-1 



0], 



(4.8) 



(4.9) 



X[t+l\t) = Mt+1) X(t\t), 
as is its eovariance matrix 



(4.1) 



the minus one appearing in the position of the clock. 
Although there are only m-1 measurements, m tests, one 
for each clock, are possible, but they are not statistically 
independent tests. Letting A represent any of the A vec- 
tors, the minimum variance weighted least squares 
estimate of () is 



P(t+\\t) = <t>(t+l) P(t\t) <J>'(t+l) + d"(t+lK>. (4.2) * = A'CT'it +l)Kt+l)/A 'C-'U+IM, (4.10) 

The predicted values of the observations at time (+1 are an< * " as standard error 

Z(t+l\t) = Hit+1) Xit+l\t). (4.3) s.e.ib) ={l/A'C-i(t+l)A) v >. (4.11) 

If the test statistic 

z = b/s.e.{b) (4.12) 



The difference between the vector of actual observa- 
tions and predicted observations is the innovation vector, 



lit+1) = Z(t+1) - Zit+l\t), (4.4) 

which has eovariance matrix 

Cit+1) = Hit+l)P(t+l]t) H'(f+1) + J?(t+1).(4.5) 

It is at this point in the recursion that a statistical 
method of error detection is employed. A simple method 
would be to divide each element of the innovation vector 
by the square root of corresponding diagonal element of 
its eovariance matrix giving a standard normal variable 
under the null hypothesis of no error. However, since clock 
differences are being measured and the elements of the 
innovation vector are intercorrelated, this is not an op- 
timal test. 

If there is an error in the reference clock, a constant 
bias will appear in every reading. This can be written as 
a regression equation with correlated errors 



has absolute value larger than some value such as 3.0, 
the corresponding clock is assumed to have an error. It 
must be remembered that this test is only approximate 
if guesses are being used for the random walk variances 
in the Q matrix, and will be much more accurate after 
maximum likelihood estimates are obtained for these 
variances. The details of the above calculations are given 
in the Appendix. 

An overall test for errors is possible by calculating the 
quadratic form 



Q = l'(t+l) Cfc+l) /(*+l| 



(4.13) 



I(t+1) = A, /J + £, 

where Aj is a column vector of ones 
AS = [1, 1, ...1], 



(4.6) 



(4.7) 



and the error vector, i, has eovariance matrix C(t+1) 
from eq (4.5). If there is an error in a clock which is not 
the reference clock, the error will appear only in the 
measurement involving that clock. For clock i, the model 



which will be distributed as chi-square with m~ 1 degrees 
of freedom under the null hypothesis of no errors in the 
observations. Although this statistic is calculated, it is no 
longer being used for error detection since it is felt that 
this test lacks power against the most common alternative 
of an error in a single clock. 

If an error is detected, the reading for the clock in ques- 
tion is eliminated from the observation vector. If more 
than one error is detected the clock with the largest ab- 
solute value of z is eliminated and the tests recalculated 
for the remaining clocks. If the clock to be eliminated is 
the reference clock, the reference clock is eliminated by 
choosing a clock still in the model to be the new reference 
clock and forming a new data vector by subtracting the 
reading of the difference between the previous reference 
clock and the new reference clock from each of the other 
readings. This is shown schematically as: 
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When all remaining clocks pass the test, the recursion 
continues. Let the Kalman gain be 

Mt+l) = P(t+l\t) ff'(H-l) CT'it+l). (4.15) 

If clocks have been eliminated due to errors, the rows 
of H{t-\-l) corresponding to the eliminated clocks are 
eliminated as are the corresponding rows and columns of 
C(t+ 1 ). In the cases where the reference clock is changed, 
the positions of the ones in the H(t+1) matrix must be 
moved to correspond to the new reference clock. The 
estimated state vector is then updated, 

X(t+l|t+l) = X(t+l\t) + A(H-l) Kt+D, (4.16) 

as is its covariance matrix 

P{t+l\t+l) = Pit+l\t) - 

Mt+1) Hit+1) Pit+l\t). (4.17) 

This completes the recursion, but the effects of any 
detected errors must be resolved. Because the measure- 
ment errors of the data are so small, the final updating 
eq (4.16) essentially sets the time states to agree with the 
actual readings of the differences. This will not be the case 
for any clock that was eliminated at this step since its 
reading has been eliminated. That clock will simply re- 
main at its position which was predicted from the previous 
time point. Since the most common type of error is a time 
step (a clock change which persists into the future), an 
administrative decision was made to correct the time state 
of a clock with a detected error so that the difference be- 
tween the time of the reference clock and the clock with 
an error agreed with the actual observation. If the error 
was a read error which does not persist into the future, 
an error will probably be detected at the next time point 
which will have opposite sign and a reverse correction 
applied. 

The remaining problem is the possibility of a frequency 
step. If a very small frequency step occurs, the Kalman 
algorithm has the flexibility to slowly adjust the estimated 
frequency of the clock to the new value. If the frequency 



step causes a time change large enough to be detected, 
it is not known in one interval whether this is a time error 
or frequency step. In case a frequency step did occur, it 
is desirable to allow the frequency state to adjust quickly 
to the new value. In order to accomplish this, the diagonal 
element of the state covariance matrix, P(t+l \t +1 ) cor- 
responding to y(r+l) of the clock in question is increased. 
The value that is added to this diagonal element is the 
minimum of (2*c/<5U)P and d(t)10 6 , where c is the clock 
time correction. This allows the recursion to adjust to a 
new frequency within a few time steps. Whether or not 
a frequency step occurred, the P matrix returns to its 
previous order of magnitude within a few time steps. 

The estimate of "true" time is obtained from the state 
vector X[t\t) which contains the estimates of x{t ) for each 
clock. The estimated values of x (t) for each clock are in 
relative agreement since their differences agree with the 
observed differences at time t. A relatively crude clock 
within the computer is used to determine the approximate 
time of the measurement. 

Each of the clocks rims freely and independently of the 
others. Thus the ticks of the clocks each second are ran- 
domly distributed. If the tick on the reference clock is used 
as the indicator to read the other clocks, the readings are 
the time differences between the ticks in nanoseconds. At 
some time in the past the ensemble was initialized relative 
to some definition of "true" time (for example, interna- 
tional Atomic Time) and the uncertainty expressed in the 
initial covariance matrix P(0|0). At that time each clock's 
tick was a certain number of nanoseconds off of "true" 
time. As data are collected and the algorithm progresses 
through its steps, all the clocks remain in relative agree- 
ment, but the ensemble drifts as a whole. The variance 
of this ensemble drift is given by the diagonal elements 
of P{t\t) corresponding to the x{t) for each clock. Since 
observational error is so small, these diagonal elements 
are essentially all the same. True time is estimated by cor- 
recting the tick of any one of the clocks by an amount 
equal to its estimated time state, x[t), which is an estimate 
of the clock's deviation from true time. The variance of 
time estimate is the corresponding diagonal element of 

Pit\t)- 

Sample output from the algorithm starting with the in- 
itialization, which was output from the previous run, is 
shown in figure 1. An example of a detected error with 
an adjustment is shown in figure 2. The algorithm also 
has the capability to adjust any clock administratively by 
a given amount or apply a frequency adjustment to all 
clocks in order to steer a time scale. These are software 
adjustments which do not affect the physical clocks but 
enable a time scale to trace another time scale such as In- 
ternational Atomic Time. The National Bureau of Stan- 
dards has several of these paper scales. 
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FIGURE 1. Example of output of algorithm starting with the initialization from the previous run.,Under INITIALIZATION and UPDATED 
STATE, the standard deviations are the square roots of the diagonal elements of the state covariance matrix, P((|f). Under MJD = 45060, 
the standard deviations are the square roots of the diagonal elements of. the innovation covariance matrix, C{t\. 
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T. 48 


0. ooo 


11 1.339 


270, a 


-96. 70 


0. 000 


a. 7 


0. 26 


o.ooo 




I4B2, 1 


B9. 33 


0.000 


14 60.1 


131B0, a 


41. 43 


0, ooo 


0. 1 


-0.02 


0. ooo 




14B2. 1 


4.62 


0. 000 


13 a 


-32123. 5 


-961 ... Bl 


0. 000 


14.8 


3. 49 


0.000 




1482. 1 


6.99 


o.ooo 



FIGURE 2. Example of computer output with' an error detected in clock 1 1 . The listing starts with the final results from the previous day. 
QUAD refers to the overall test which is not actually used for detecting errors but is printed out so that errors are easily visible. The optimal 
estimate of. the clock error is -44.7 with standard deviation 11.7. 
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5. Estimation of Parameters 

Assuming Gaussian errors, -2 In likelihood is calculated 
from 



CLOCKS PRESENT IN STUDY 



L = Z[/n|C(t)| + Tit) C-'dl/W]. 
t 



(5.1) 



This gives a number which depends on the values of 
the parameters of the model. Nonlinear optimization 
routines can then be used to minimize this function giv- 
ing maximum likelihood estimates of the parameters. 
Each function evaluation is a pass of the recursion through 
all the data. 

When error detection is used, a new problem arises. The 
dimensions of C(t) and I{t) are reduced if a clock is 
eliminated at time t. When the nonlinear optimization 
routine varies the parameters, it is possible that at some 
point a clock that was previously eliminated is suddenly 
not eliminated, or that a clock that was not eliminated 
is eliminated. This would cause a discontinuity in the 
likelihood function. To avoid this problem, a file of in- 
dicator variables, indicating which clocks are deleted at 
which time points, is produced. The nonlinear optimiza- 
tion is then run conditional on fixed clocks being 
eliminated at certain time points. Upon convergence the 
parameters in the original model are replaced by their op- 
timized value and the initial program rerun generating a 
new array of indicator variables. If the vaiue of the 
likelihood function changes it indicates that the error 
detection procedure has detected a different set of errors 
so the nonlinear optimization is repeated. Experience has 
shown this procedure converges in two or three iterations. 

The optimization was carried out on one year of data, 
April 1, 1981 to March 31, 1982. During this time 12 
clocks were used, but since clocks were added and deleted 
administratively from the algorithm (removed because of 
problems with the clock, not error detection) a maximum 
of 10 clocks and a minimum of eight clocks were available 
at any one time. Figure 3 shows the times at which various 
clocks were present and the number of actual observa- 
tions used for each clock. The number missing is the 
number of observations eliminated by error detection. 

In the first nonlinear optimization run, o £ and o„ were 
varied for each clock and w(t) and o a were set equal to 
zero. This is a drift-free model, and with 12 clocks gives 
24 parameters to be estimated. Standard deviations rather 
than variances are used for optimization since if they go 
negative their squares are still positive. If variances were 
used directly, it would be possible for the nonlinear op- 
timization routine to try negative values which would be 
meaningless. The results of this optimization procedure 
are shown in table 1. 



Clack 
601 

B — 
113 — 

NBS4 — 

1375 — 
323 — 
352 — 

PHM4 

SI — 

137 — 

167 — 

1316 — 




No. 

Pi-WBimt 

29B 


Nd. 

Mil. ulna 

2 


360 


5 


354 


5 


56 


3 


357 


B 


255 


5 


354 


11 


203 


6 


67 


4 


35B 


7 


361 


4 


364 


1 



365 



FIGURE 3. A schematic diagram showing the time when various clocks 
are available for use in the algorithm. The number present are the actual 
number used after error detection. The numbers missing are the number 
of errors detected. The sum of the number present and number missing 
is the number of observations available to the algorithm. 



Table 1 . Estimated values of o, and o, an d 95% confidence intervals. 



Clock 


n 


a. 


"j nanoseconds 


o, ™ 


nanoseconds/day 






Lower 




Upper 


Lower 




Upper 






Limit 


Est 


Limit 


Limit 


Est 


Limit 


1316 


364 


3.81 


4.14 


4.53 


0.53 


0.80 


1.23 


167 


361 


12.58 


13.52 


14.67 


0.57 


1.11 


2.07 


137 


358 


10.41 


11.31 


12.27 


1.76 


2.49 


3.56 


61 


67 


5.48 


6.77 


8.43 


1.53 


2.80 


4.83 


352 


354 


8.12 


8.85 


9.74 


2.42 


3.32 


4.41 


323 


255 


2.06 


2,37 


2.74 


0.63 


0.94 


1.34 


1375 


357 


9.93 


10.71 


11.64 


0.96 


1.48 


2.25 


NBS4 


66 





0.88 


1.86 


0.72 


1.34 


2.16 


113 


354 


8.73 


9.48 


10.38 


2.49 


3.18 


4.11 


8 


360 


7.98 


8.65 


9.49 


2.11 


2.76 


3.66 


601 


298 


1.89 


2.13 


2.41 





0.06 


0.52 


PHM4 


203 





0.65 


1.19 


0.55 


0.77 


1.09 



When m clocks are in the model, it is only possible to 
estimate m-1 drift terms when the observations are time 
differences. Tryon and Jones [1] handled this problem by 
estimating m drifts with the constraint that the sum of 
the drifts was zero. Here, a different approach is used. 
A clock in the model is assumed to have zero drift and 
the drifts of the m-1 remaining drifts are estimated. Clock 
601 was chosen to be the clock with no drift since the stan- 
dard deviation of its frequency random walk was not 
significantly different from zero. Since clock 601 was not 
available for the entire year, the drift for clock 323 was 
also set equal to zero. These two clocks had overlapping 
data intervals which included the entire year. Clock 323 
was also chosen because it is known to be a particularly 
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stable clock. The drifts for two other clocks, 61 and NBS4, 
were also set equal to zero because of their short lengths 
of data. The remaining eight clocks were tested for drift 
by a 24-parameter nonlinear optimization run. For each 
of the eight clocks, the frequency random walk standard 
deviation, o„, the initial value of the drift, w{0), and the 
drift random walk standard deviation, a a , were allowed 
to vary. All of the time random walk standard deviations, 
o £ , were fixed at their optimized values from the previous 
run, as were the parameters for the remaining four clocks 
in the algorithm. This approach of fixing some of the 
parameters was necessitated by the storage available on 
the PDP 11/70 which allowed a maximum of 25 
parameters to be optimized simultaneously. Changes in 
o £ that occur when drift is included in the model, while 
small, do in fact exist [1]. 

The estimated values of all eight o a were very small, 
so their significance was tested by fixing their values at 
zero and rerunning the previous optimization with 16 
parameters. Statistical tests based on changes in -2 In 
likelihood have asymptotic chi-square distributions with 
degrees of freedom equal to the number of parameters 
added or deleted from the model or constraints introduced 
[2]. In this case, eight parameters were fixed at zero but 
-2 In likelihood only increased by 0.2 indicating that the 
a a are not significantly different from zero. This agrees 
with the results obtained on the earlier data set [1], These 
results, which were obtained on a CDC Cyber 170/750 
computer, were checked using the new algorithm on the 
PDP 11/70 and duplicated exactly. At this point it is 
assumed that any frequency drifts which may exist can 
be considered to be constant over the length of the data 
span. 

When the eight drift terms were eliminated from the 
model, the change in -2 Ire likelihood was not significant, 
indicating that the inclusion of drift terms did not improve 
the model significantly. This did not agree with the 
previous study where there were significant trends for 
clocks 601, 137, and 61, but there are explanations for 
this. In the first study clock 601 was nearing the end of 
its life, and the cesium beam tube was replaced between 
the two studies making it essentially a different clock. 
Clock 137 was introduced as a new clock on the first day 
of the first study. Clocks are known to drift near the begin- 
ning and end of their lives and be more stable during 
midlife. Clock 61 was only available for 69 days in the 
second study so any drift that may have existed could not 
be detected. If a common drift exists in all clocks, this 
could not be estimated, and the uncertainty introduced 
into the time scale would not appear in the P(t|t) matrix. 

It is of interest to compare the confidence intervals 
obtained on two different sets of data. Table 2 shows the 



results for the clocks that were common to the two studies. 
The estimates and confidence intervals are in agreement 
except for clock 601 in which the cesium beam tube was 
replaced. 



Table 2. 1919 and 1981 clock characterization 
95% confidence intervals. 



Clock Study o. 


^ nanoseconds 


c, "-■ 


nanoseconds/day 






Lower 




Upper 


Lower 




Upper 






Limit 


Est 


Limit 


Limit 


Est 


Limit 


601 


1 


6.88 


7.46 


8.14 





0.44 


0.96 




2 


1.89 


2.13 


2.41 





0.06 


0.52 


167 


1 


12.44 


13.45 


14.61 


0.56 


1.11 


1.96 




2 


12.58 


13.52 


14.67 


0.57 


1.11 


2.07 


137 


1 


9.26 


10.04 


10.96 


1.03 


1.60 


2.26 




2 


10.41 


11.31 


12,27 


1.76 


2.49 


3.56 


1316 


1 


3.17 


3.62 


4.05 


0.97 


1.36 


1.82 




2 


3.81 


4.14 


4.53 


0.53 


0,80 


1.23 


323 


1 


3.11 


3.53 


3.94 


0.41 


0.73 


1.10 




2 


2.06 


2.37 


2.74 


0.63 


0.94 


1.34 


8 


1 


8.32 


9.09 


9.95 


2.01 


2.65 


3.48 




2 


7.98 


8.65 


9.49 


2.11 


2.76 


3.66 



6. Conclusion 

A major advantage of a Kalman algorithm for 
estimating time is its flexibility when changes are needed. 
An example is the introduction of frequency calibrations 
from a frequency standard. A frequency calibration, 
together with its precision, can easily be introduced into 
the state vector and state covariance matrix. Experimen- 
tation with different Kalman models are possible since 
changing the state transition matrix is not difficult. 
However, the major advantage is that the Kalman 
algorithm is based on a mathematical model which fits 
the data. The innovations can be tested for whiteness, and 
standard statistical methods can be used for obtaining con- 
fidence intervals and hypothesis tests. 
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Appendix 

For computational efficiency and numerical stability, 
a procedure involving a Cholesky decomposition is used 
in the Kalman recursion [4]. The tests for errors are also 
simplified. The matrix product HU+1) Pit+l\t) in eq. 
(4.5) is calculated without actually performing a matrix 
multiplication. Since each row of Hit+l) contains a one 
and a minus one with the rest of the elements zero, each 
element of the matrix product is simply the difference be- 
ween two elements of Pit+l\t). The matrix C(t+l) is 
calculated and augmented as follows 



lCit+l):fflt+l)Pit+l\t\\Iit+-l) : .A 1 iA2: ... :A n 
(m-1) 3m 111 — 1 



(A.l) 



This matrix has m-I rows with the number of columns 
shown under each partition. Only the upper triangular 
position of dt+1) need be stored. The third partition is 
the vector of innovations from eq (4.4) followed by the 
vectors to be used for testing for clock errors. One of these 
vectors will correspond to the reference clock eq (4.7) and 
the others to the remaining clocks eq (4.9). 



The Cholesky decomposition factors C(t+1) into 

Cit+1) = T'it+1) Tit+l), (A.2) 

where Tit + 1 ) is upper triangular. The algorithm as given 
in [4] works on the entire augmented matrix and is 
equivalent to premultiplying the matrix by [T'(H-1)]~ L . 
Let 

Bit+1) = [T'it+U]- 1 Hit+l) Pit+l\t) 
DiT+l) = IT'U+l)]- 1 Iit+1) (A.3) 

G t = [T'U+l]]- 1 A { . 
The test for each clock, eq (4.10) becomes 

b ■ = G'Dit+D/G'G, (A.4) 

and eq (4.11) becomes 

8jn.[b) = (l/G'G,w (A.5) 

The quadratic form in eq (4.13) reduces to 

Q = D'it+l) DU+1), (A.6) 

and eq (4.15) to 

A(H-l) = B'it+1) [T'it+DY 1 . (A.7) 

The final updating equations are (4.16) 
XU+l\t+l) = X(r+l|r) 4- B'[t+1) D(t+1), (A.8) 
and (4.17) 
P(t+l|t+l) = J>U+l|t) - B'(t+1) B(t+1). (A.9) 



24 



